Nonlinear mixed-effect models for individual differences in growth

Load libraries

library(tidyverse); library(patchwork); library(tidybayes)
library(brms); library(marginaleffects); library(viridis)
library(easystats)

Set theme

theme_set(theme_bw(16))

Model importation

# Importing all models from \stan folder to avoid refitting
# Note, there must be a way to automatically import all models files
# from a folder. Doing it the hard way for now
gomp.prior.default = readRDS("stan/gomp.prior.default.rds")
gomp.prior = readRDS("stan/gomp.prior.rds")
gomp = readRDS("stan/gomp.rds")
gomp.sc.prior = readRDS("stan/gomp.sc.prior.rds")
gomp.sc = readRDS("stan/gomp.sc.rds")
gomp.to = readRDS("stan/gomp.to.rds")

Rationale

I’ve been working a lot with nonlinear models lately. In many cases across ecology and evolution, trends do not go on forever. Having modeling approaches that can account for boundaries and plateaus is therefore very important in specific situations. In my case I’m generally interested in how contamination exposure affects growth, reproduction and behavior of invertebrates and its consequences on individual fitness. In this post I’m exploring how to use nonlinear mixed effect models of growth to recover individual trajectories using simulated data. I’m mostly using the {brms} package for now as it provides a very user-friendly syntax for fitting these models. The level of complexity I’m dealing with is veering dangerously toward a full STAN implementation, especially when I’ll have to deal with jointly modeling fitness, but that’s a problem for future me!

Inspiration

A lot of the code in this document has been taken and adapted from Solomon Kurtz’s work and his translation of Richard McElreath’s Statistical Rethinking book into {brms} syntax. Specifically, I’ve been working through Chapter 14 section 1.3.Many thanks to them for putting out all this material in an open format!

Model description

I’m using the model described in Goussen et al. (2013) showing adaptation of C. elegans to Uranium. The model assumes nematodes body-length grows according to the following Gompertz function:

\[ L = L_{inf}\times e^{ln \left( \frac{L_0}{Linf} \right) \times e^{-\alpha t} }\]

where \(L_{inf}\) is the asymptotic body-length, \(L_0\) is the body-length at hatching and \(\alpha\) the growth rate.

Plugging in the values from the paper gives a sigmoid growth curve:

# Store parameters for non-exposed populations
n_id = 10 # 10 individuals
obs = seq(0, 126, by = 24) # One observation every day
Linf = 1370
L0 = 182
alpha = 0.028
t = seq(0,126, by = 1) # measure every h for 126 h 
sigma = 100 # random noise

Gomp.fun = function(t, Linf, L0, alpha){
  Lhat = Linf * exp(log(L0/Linf)*exp(-alpha*t)) # predicted growth
  return(Lhat)
}

# Apply equation
set.seed(42)
df = crossing(obs, 1:n_id) %>% 
  mutate(Lhat = Gomp.fun(obs, Linf = Linf, L0 = L0, alpha = alpha)) %>% 
  mutate(L = rnorm(n(), Lhat, sigma))
df %>% 
  ggplot(aes(y = L, x = obs)) +
  geom_point(alpha = .2, size = 2.5) +
  geom_function(fun = Gomp.fun, 
                     args = list(Linf = Linf, L0 = L0, alpha = alpha),
                     color = "red", size = 1) +
  ylim(0, 1600) +
  labs(x = "Hours since hatching", y = "Body-length (mum)")

Simulating individual differences

We can use the brms package to make simulations according to some prior distributions for our parameters. We first need to store the values of the growth parameters for each individual

Data simulations

I’m using the same approach as above but now allow individual to have some deviation from the population mean for all parameters: \(\mu_i \sim N(mu + \text{offset}, \text{ } \sigma)\) where \(\mu\) and \(\mu_i\) are the individual and population mean for a given parameter and the offset is calculated as: \(\text{offset} \sim N(0, \text{ } \sigma_i)\) where \(\sigma_i\) is the among-individual standard deviation for a given parameters. In our simulation, we assume a coefficient of variation 10 % for all parameters such that \(\sigma_i = 0.10 \times \mu\)

Note

There are different ways to interpret the \(L_{inf}\) parameter. In Goussen et al’ paper, this was assumed to be the maximum body-size for the population. Meaning that with enough time, all individuals will converge toward this value. Most individuals of course will die before that! Here, I am assuming that all individuals have their own maximum size they can reach. This means that two individuals with the same growth rate measured for the same amount of time may still differ in their maximum length.

n_id = 30 # 30 individuals
times = seq(0, 126, by = 24) # One observation every day
sd = 10 # random noise
L0_mu = 182 # initial length (micrometers)
Linf_mu = 1370 # maximal length (micrometers)
alpha_mu = 0.028 # Growth rate (hour^-1)

rho = 0 # Suppose all parameters are independent

mu     = c(L0_mu, Linf_mu, alpha_mu)
sigmas = c(L0_mu*.1, Linf_mu*.1, alpha_mu*.1) # 10 % CV around the mean
rho_mat = matrix(c(1, rho, rho,
                    rho, 1, rho,
                    rho, rho, 1), 
                    nrow = 3)

sigma = diag(sigmas) %*% rho_mat %*% diag(sigmas)

set.seed(42)
ID = MASS::mvrnorm(n_id, mu, sigma) %>% 
  data.frame() %>% 
  set_names("L0_i", "Linf_i", "alpha_i") %>% 
  mutate(ID = 1:n_id)

# Simulate individual growth
df = ID %>%
  tidyr::expand(nesting(ID, L0_i, Linf_i, alpha_i), 
                t = times) %>%
  mutate(Lhat = Linf_i * exp(log(L0_i/Linf_i)*exp(-alpha_i*t))) %>%
  mutate(L = rnorm(n(), Lhat, sd))

df %>% 
  ggplot(aes(y = L, x = t, group = ID)) +
  geom_point(alpha = .2, size = 2.5) +
  geom_line(alpha = .2, size = .5) +
  geom_function(fun = Gomp.fun, 
                     args = list(Linf = Linf, L0 = L0, alpha = alpha),
                     color = "red", size = 1) +
  ylim(0, 1600) +
  labs(x = "Hours since hatching", y = "Body-length (mum)")

Defining the model in brms

In brms, we can use the standard R formula syntax to specify the Gompertz model

bf = bf(L ~ 
           Linf * exp(log(L0 / Linf) * exp(-alpha * t)), # Gompertz population curve
         L0 + Linf + alpha ~ 1 + (1|ID), # parameters with random effects
         nl = T)
priors = get_prior(bf, df)
priors

                  prior class      coef group resp dpar nlpar lb ub
 student_t(3, 0, 521.9) sigma                                  0   
                 (flat)     b                           alpha      
                 (flat)     b Intercept                 alpha      
 student_t(3, 0, 521.9)    sd                           alpha  0   
 student_t(3, 0, 521.9)    sd              ID           alpha  0   
 student_t(3, 0, 521.9)    sd Intercept    ID           alpha  0   
                 (flat)     b                              L0      
                 (flat)     b Intercept                    L0      
 student_t(3, 0, 521.9)    sd                              L0  0   
 student_t(3, 0, 521.9)    sd              ID              L0  0   
 student_t(3, 0, 521.9)    sd Intercept    ID              L0  0   
                 (flat)     b                            Linf      
                 (flat)     b Intercept                  Linf      
 student_t(3, 0, 521.9)    sd                            Linf  0   
 student_t(3, 0, 521.9)    sd              ID            Linf  0   
 student_t(3, 0, 521.9)    sd Intercept    ID            Linf  0   
       source
      default
      default
 (vectorized)
      default
 (vectorized)
 (vectorized)
      default
 (vectorized)
      default
 (vectorized)
 (vectorized)
      default
 (vectorized)
      default
 (vectorized)
 (vectorized)

These are default priors and may not be well-suited for our problem. For some reason, brms won’t let us sample from these default priors, but we can easily illustrate the point by setting unrealistically uninformative priors (i.e. normal(0, 100) here) on these parameters. I’m using the lower bound argument lb = 0 to keep the growth positive.

priors = priors = 
  # Intercept priors
  prior(normal(0, 100), nlpar = L0, class = b, lb = 0) +
  prior(normal(0, 100), nlpar = Linf, class = b, lb = 0) +
  prior(normal(0, 100), nlpar = alpha, class = b, lb = 0)

Next, we fit the model and simply specify sample_prior = "only" in the brm() function to only get the growth trends implied by the priors

# plotting
plot(conditional_effects(gomp.prior.default, ndraws = 100, spaghetti = T))

Not very good right?! Let’s retry now with more reasonable priors

priors = 
  # Intercept priors
  prior(normal(200, 40), nlpar = L0, class = b, lb = 0) +
  prior(normal(1500, 300), nlpar = Linf, class = b, lb = 0) +
  prior(normal(.03,.006), nlpar = alpha, class = b, lb = 0) + 
  # Random effects priors (informative priors with 20 % CV)
  prior(exponential(.025), nlpar = L0, class = sd, group = ID) +
  prior(exponential(.003), nlpar = Linf, class = sd, group = ID) +
  prior(exponential(170), nlpar = alpha, class = sd, group = ID) +
  # Residual prior
  prior(exponential(1), class = sigma) 

# Plot priors
p1 = priors %>% 
  parse_dist() %>% 
  filter(class == "b") %>% 
  ggplot(aes(xdist = .dist_obj, y = format(.dist_obj))) +
  stat_dist_halfeye() +
  facet_wrap(~nlpar, scales = "free") +
  ggtitle("Intercepts") +
  xlab("Value") + ylab("Density") +
  theme_bw(12) +
  theme(axis.text.y = element_text(angle = 90)) 

p2 = priors %>% 
  parse_dist() %>% 
  filter(class == "sd") %>% 
  ggplot(aes(xdist = .dist_obj, y = format(.dist_obj))) +
  stat_dist_halfeye() +
  facet_wrap(~nlpar, scales = "free") +
  ggtitle("Among-individual variance") +
  xlab("Value") + ylab("Density") +
  theme_bw(12) +
  theme(axis.text.y = element_text(angle = 90)) 

p3 = priors %>% 
  parse_dist() %>% 
  filter(class == "sigma") %>% 
  ggplot(aes(xdist = .dist_obj, y = format(.dist_obj))) +
  stat_dist_halfeye() +
  facet_wrap(~nlpar, scales = "free") +
  ggtitle("Residual variance") +
  xlab("Value") + ylab("Density") +
  theme_bw(12) +
  theme(axis.text.y = element_text(angle = 90)) 
 
(p1 + p2 + p3) + plot_layout(ncol = 1)

Prior predictive checks

We can now check if the model shows a more appropriate gortwh curve while sampling only from these new priors.

Inspecting the model and plotting predicted growth

gomp.prior

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: L ~ Linf * exp(log(L0/Linf) * exp(-alpha * t)) 
         L0 ~ 1 + (1 | ID)
         Linf ~ 1 + (1 | ID)
         alpha ~ 1 + (1 | ID)
   Data: df (Number of observations: 180) 
  Draws: 4 chains, each with iter = 1000; warmup = 500; thin = 1;
         total post-warmup draws = 2000

Group-Level Effects: 
~ID (Number of levels: 30) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(L0_Intercept)       40.67     42.23     1.14   156.94 1.00     2507     1141
sd(Linf_Intercept)    332.87    334.35     8.69  1223.07 1.00     2706      955
sd(alpha_Intercept)     0.01      0.01     0.00     0.02 1.00     2892     1169

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
L0_Intercept      198.73     41.39   120.14   281.26 1.00     3330     1271
Linf_Intercept   1492.43    311.92   885.14  2091.52 1.00     3261     1357
alpha_Intercept     0.03      0.01     0.02     0.04 1.00     3684     1428

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.03      1.00     0.04     3.46 1.00     2297     1327

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(conditional_effects(gomp.prior, ndraws = 100, spaghetti = T))

re = crossing(t = seq(min(df$t), 
                       max(df$t),
                       length.out=100),
               ID = unique(df$ID)) %>% 
  add_epred_draws(gomp.prior, re_formula = NULL, scale = "response", ndraws = 20)

re %>% 
  ggplot(aes(y = .epred, x = t)) +
  geom_line(aes(y = .epred, x = t, group = .draw), size = .5, alpha = .5) +
  geom_point(data = df, aes(y=L, x=t, color = ID)) +
  facet_wrap(~ID, nrow = 6, ncol = 5) + 
  scale_color_viridis() +
  ylab("Mass (mg)") + 
  xlab("Time") +
  theme_bw(12) +
  theme(legend.position = "none")

Fit model to simulated data

Next, we can fit the model to the simulated dataset

Inspecting the model and plotting predicted growth

model_parameters(gomp, effects = "all")

# Fixed Effects

Parameter       |  Median |            95% CI |   pd |   Rhat |  ESS
--------------------------------------------------------------------
L0_Intercept    |  179.45 | [170.68,  453.99] | 100% | 16.769 | 2.00
Linf_Intercept  | 1362.79 | [  0.08, 1434.54] | 100% | 23.184 | 2.00
alpha_Intercept |    0.03 | [  0.00,    0.03] | 100% | 32.253 | 2.00

# Sigma

Parameter | Median |         95% CI |   pd |   Rhat |  ESS
----------------------------------------------------------
sigma     |   9.75 | [8.36, 139.79] | 100% | 20.582 | 2.00

# Random Effects Variances

Parameter                |   Median |         95% CI |   pd |  Rhat |  ESS
--------------------------------------------------------------------------
SD (L0_Intercept: ID)    |    18.31 | [1.50,  27.59] | 100% | 1.228 | 6.00
SD (Linf_Intercept: ID)  |   156.11 | [0.01, 215.75] | 100% | 3.979 | 2.00
SD (alpha_Intercept: ID) | 2.24e-03 | [0.00,   0.00] | 100% | 1.515 | 3.00

pp_check(gomp, ndraws = 100)

plot(conditional_effects(gomp, ndraws = 100, spaghetti = T))

re = crossing(t = seq(min(df$t), 
                       max(df$t),
                       length.out=100),
               ID = unique(df$ID)) %>% 
  add_epred_draws(gomp, re_formula = NULL, scale = "response", ndraws = 20)

re %>% 
  ggplot(aes(y = .epred, x = t)) +
  geom_line(aes(y = .epred, x = t, group = .draw), size = .5, alpha = .5) +
  geom_point(data = df, aes(y=L, x=t, color = ID)) +
  facet_wrap(~ID, nrow = 6, ncol = 5) + 
  scale_color_viridis() +
  ylab("Mass (mg)") + 
  xlab("Time") +
  theme_bw(12) +
  theme(legend.position = "none")

The model functions correctly but there are a lot of divergent transitions and some poor convergence diagnostics. Let’s try reparametrizing to improve the sampling.

Model reparametrization

We can express all body-length parameters as a function of the maximum length to help the model converge. Divinding each side of the Gompertz equation gives

\[\begin{aligned} l &= L/L_{inf}\\ l_0 &=L_0/L_{inf}\\ l &= e^{ln \left( l_0 \right) \times e^{-\alpha t} } \end{aligned}\]

ggplot() +
  geom_function(fun = Gomp.fun, 
                     args = list(Linf = 1, L0 = L0/Linf, alpha = alpha),
                     color = "red", size = 1) +
  xlim(0, 126) + ylim(0, 1) +
  labs(x = "Hours since hatching", y = "Scaled body-length") 

Prior predictive checks

We need to reformulate model formula and priors in order to get sensible outputs

df$l = with(df, L/Linf)

bf = bf(l ~ 
           exp(log(l0) * exp(-alpha * t)), # Gompertz population curve
         l0 + alpha ~ 1 + (1|ID), # parameters with random effects
         nl = T)

Defining and plotting priors

priors = 
  # Intercept priors
  prior(normal(.15, .03), nlpar = l0, class = b, lb = 0) +
  prior(normal(.03,.006), nlpar = alpha, class = b, lb = 0) + 
  # Random effects priors (informative priors with 20 % CV)
  prior(exponential(34), nlpar = l0, class = sd, group = ID) +
  prior(exponential(170), nlpar = alpha, class = sd, group = ID) +
  # Residual prior
  prior(exponential(1), class = sigma) 

# Plot priors
p1 = priors %>% 
  parse_dist() %>% 
  filter(class == "b") %>% 
  ggplot(aes(xdist = .dist_obj, y = format(.dist_obj))) +
  stat_dist_halfeye() +
  facet_wrap(~nlpar, scales = "free") +
  ggtitle("Intercepts") +
  xlab("Value") + ylab("Density") +
  theme_bw(12) +
  theme(axis.text.y = element_text(angle = 90)) 

p2 = priors %>% 
  parse_dist() %>% 
  filter(class == "sd") %>% 
  ggplot(aes(xdist = .dist_obj, y = format(.dist_obj))) +
  stat_dist_halfeye() +
  facet_wrap(~nlpar, scales = "free") +
  ggtitle("Among-individual variance") +
  xlab("Value") + ylab("Density") +
  theme_bw(12) +
  theme(axis.text.y = element_text(angle = 90)) 

p3 = priors %>% 
  parse_dist() %>% 
  filter(class == "sigma") %>% 
  ggplot(aes(xdist = .dist_obj, y = format(.dist_obj))) +
  stat_dist_halfeye() +
  facet_wrap(~nlpar, scales = "free") +
  ggtitle("Residual variance") +
  xlab("Value") + ylab("Density") +
  theme_bw(12) +
  theme(axis.text.y = element_text(angle = 90)) 
 
(p1 + p2 + p3) + plot_layout(ncol = 1)

And now running the model on priors only

Model inspection

model_parameters(gomp.sc.prior, effects = "all")

# Fixed Effects

Parameter       | Median |       95% CI |   pd |  Rhat |     ESS
----------------------------------------------------------------
l0_Intercept    |   0.15 | [0.09, 0.20] | 100% | 1.001 | 2770.00
alpha_Intercept |   0.03 | [0.02, 0.04] | 100% | 1.000 | 2270.00

# Sigma

Parameter | Median |       95% CI |   pd |  Rhat |     ESS
----------------------------------------------------------
sigma     |   0.69 | [0.02, 3.91] | 100% | 1.000 | 2914.00

# Random Effects Variances

Parameter                |   Median |       95% CI |   pd |  Rhat |     ESS
---------------------------------------------------------------------------
SD (l0_Intercept: ID)    |     0.02 | [0.00, 0.12] | 100% | 0.998 | 2925.00
SD (alpha_Intercept: ID) | 4.03e-03 | [0.00, 0.02] | 100% | 0.999 | 2851.00

plot(conditional_effects(gomp.sc.prior, ndraws = 100, spaghetti = T))

re = crossing(t = seq(min(df$t), 
                       max(df$t),
                       length.out=100),
               ID = unique(df$ID)) %>% 
  add_epred_draws(gomp.sc.prior, re_formula = NULL, scale = "response", ndraws = 20)

re %>% 
  ggplot(aes(y = .epred, x = t)) +
  geom_line(aes(y = .epred, x = t, group = .draw), size = .5, alpha = .5) +
  geom_point(data = df, aes(y=l, x=t, color = ID)) +
  facet_wrap(~ID, nrow = 6, ncol = 5) + 
  scale_color_viridis() +
  ylab("Mass (mg)") + 
  xlab("Time") +
  theme_bw(12) +
  theme(legend.position = "none")

Fitting model to data

As above we rerun the scaled model on the scaled data

Model inspection

model_parameters(gomp.sc, effects = "all")

# Fixed Effects

Parameter       | Median |       95% CI |   pd |  Rhat |     ESS
----------------------------------------------------------------
l0_Intercept    |   0.13 | [0.12, 0.14] | 100% | 1.000 | 4283.00
alpha_Intercept |   0.03 | [0.03, 0.03] | 100% | 1.004 |  658.00

# Sigma

Parameter | Median |       95% CI |   pd |  Rhat |     ESS
----------------------------------------------------------
sigma     |   0.04 | [0.03, 0.04] | 100% | 1.000 | 3774.00

# Random Effects Variances

Parameter                |   Median |       95% CI |   pd |  Rhat |    ESS
--------------------------------------------------------------------------
SD (l0_Intercept: ID)    | 6.23e-03 | [0.00, 0.02] | 100% | 1.008 | 989.00
SD (alpha_Intercept: ID) | 5.63e-03 | [0.00, 0.01] | 100% | 1.003 | 733.00

plot(conditional_effects(gomp.sc, ndraws = 100, spaghetti = T))

pp_check(gomp.sc, ndraws = 100)

re = crossing(t = seq(min(df$t), 
                       max(df$t),
                       length.out=100),
               ID = unique(df$ID)) %>% 
  add_epred_draws(gomp.sc, re_formula = NULL, scale = "response", ndraws = 20)

re %>% 
  ggplot(aes(y = .epred, x = t)) +
  geom_line(aes(y = .epred, x = t, group = .draw), size = .5, alpha = .5) +
  geom_point(data = df, aes(y=l, x=t, color = ID)) +
  facet_wrap(~ID, nrow = 6, ncol = 5) + 
  scale_color_viridis() +
  ylab("Mass (mg)") + 
  xlab("Time") +
  theme_bw(12) +
  theme(legend.position = "none")

Much better!

Assuming a trade-off between maximum size and growth-rate

So far, we have assumed that all our parameters are independent. However, trade-offs between various life-history traits are also likely to occur. One such trade-off is the compromise between fast growth and reproductive rate, such that individuals that mature sooner reach a smaller asymptotic length. This is the case in collembolas exposed to Cadmium for example. We can easily modify our simulation function to incorporate a similar trade-off in our data.

Data simulations and modeling

n_id = 30 # 30 individuals
times = seq(0, 126, by = 24) # One observation every day
sd = 10 # random noise
L0_mu = 182 # initial length (micrometers)
Linf_mu = 1370 # maximal length (micrometers)
alpha_mu = 0.028 # Growth rate (hour^-1)

rho = -.7 # Assume strong negative correlation between Linf and alpha

mu     = c(L0_mu, Linf_mu, alpha_mu)
sigmas = c(L0_mu*.1, Linf_mu*.1, alpha_mu*.1) # 10 % CV around the mean
rho_mat = matrix(c(1, 0, 0,
                    0, 1, rho,
                    0, rho, 1), 
                    nrow = 3)

sigma = diag(sigmas) %*% rho_mat %*% diag(sigmas)

set.seed(42)
ID = MASS::mvrnorm(n_id, mu, sigma) %>% 
  data.frame() %>% 
  set_names("L0_i", "Linf_i", "alpha_i") %>% 
  mutate(ID = 1:n_id)

# Plot
ID %>% 
  select(-ID) %>% 
  GGally::ggpairs() +
  theme_bw() 

And now we simulate growth over 126 hours as before

# Simulate individual growth
df = ID %>%
  tidyr::expand(nesting(ID, L0_i, Linf_i, alpha_i), 
                t = times) %>%
  mutate(Lhat = Linf_i * exp(log(L0_i/Linf_i)*exp(-alpha_i*t))) %>%
  mutate(L = rnorm(n(), Lhat, sd))

df %>% 
  ggplot(aes(y = L, x = t, group = ID)) +
  geom_point(alpha = .2, size = 2.5) +
  geom_line(alpha = .2, size = .5) +
  geom_function(fun = Gomp.fun, 
                     args = list(Linf = Linf_mu, L0 = L0_mu, alpha = alpha_mu),
                     color = "red", size = 1) +
  ylim(0, 1600) +
  labs(x = "Hours since hatching", y = "Body-length (mum)")

Switching directly to fitting the model to the data, we make some small adjustments to the formula

df$l = with(df, L/Linf_mu)

bf = bf(l ~ 
           linf*exp(log(l0/linf) * exp(-alpha * t)), # Gompertz population curve
         l0 + linf + alpha ~ 1 + (1|c|ID), # parameters with random effects
         nl = T)

We introduce a scaled version of \(L_{inf}\) into the formula and allow random effects to be correlated with the (1|c|ID) chunk. Next we need to specify a strong prior for on \(l_{inf}\) to constrain it around 1. We also specify a prior on the correlation parameter using an LKJ prior of parameter \(\eta>1\)

priors = 
  # Intercept priors
  prior(normal(.15, .03), nlpar = l0, class = b, lb = 0) +
  prior(normal(1, .03), nlpar = linf, class = b, lb = 0) +
  prior(normal(.03,.006), nlpar = alpha, class = b, lb = 0) + 
  # Random effects priors (informative priors with 20 % CV)
  prior(exponential(34), nlpar = l0, class = sd, group = ID) +
  prior(exponential(5), nlpar = linf, class = sd, group = ID) +
  prior(exponential(170), nlpar = alpha, class = sd, group = ID) +
  # Residual prior
  prior(exponential(1), class = sigma) +
  # Correlation prior
  prior(lkj(2), class = cor)


# Plot priors for linf
priors %>% 
  parse_dist() %>% 
  filter(nlpar == "linf") %>% 
  ggplot(aes(xdist = .dist_obj, y = format(.dist_obj))) +
  stat_dist_halfeye() +
  facet_wrap(~nlpar, scales = "free") +
  ggtitle("Scaled asymptotic length priors") +
  xlab("Value") + ylab("Density") +
  theme_bw(12) +
  theme(axis.text.y = element_text(angle = 90)) 

To get a better sense of what these prior imply, we can simulate the distribution of individual differences in \(l_{inf}\)

linf_sim = data.frame(offsets = rnorm(n_id*100, 0, .2)) %>% 
  mutate(linf_mu = rnorm(n(), (1 + offsets), sigma))
  
linf_sim %>% 
  ggplot(aes(x = linf_mu)) +
  stat_halfeye() +
  xlim(-2,3)

We now fit the model to the simulated data

Model inspection

model_parameters(gomp.to, effects = "all")

# Fixed Effects

Parameter       | Median |       95% CI |   pd |  Rhat |     ESS
----------------------------------------------------------------
l0_Intercept    |   0.13 | [0.12, 0.14] | 100% | 1.000 | 2221.00
linf_Intercept  |   1.00 | [0.97, 1.04] | 100% | 1.001 | 1854.00
alpha_Intercept |   0.03 | [0.03, 0.03] | 100% | 1.000 | 1866.00

# Sigma

Parameter |   Median |       95% CI |   pd |  Rhat |     ESS
------------------------------------------------------------
sigma     | 7.16e-03 | [0.01, 0.01] | 100% | 1.000 | 1479.00

# Random Effects Variances

Parameter                                |   Median |         95% CI |     pd |  Rhat |     ESS
-----------------------------------------------------------------------------------------------
SD (l0_Intercept: ID)                    |     0.01 | [ 0.01,  0.02] |   100% | 1.000 | 2074.00
SD (linf_Intercept: ID)                  |     0.12 | [ 0.09,  0.15] |   100% | 1.001 | 1507.00
SD (alpha_Intercept: ID)                 | 2.84e-03 | [ 0.00,  0.00] |   100% | 1.002 | 1641.00
Cor (l0_Intercept~linf_Intercept: ID)    |    -0.33 | [-0.62,  0.02] | 96.55% | 1.004 | 1004.00
Cor (l0_Intercept~alpha_Intercept: ID)   |     0.07 | [-0.30,  0.43] | 64.42% | 1.001 | 1036.00
Cor (linf_Intercept~alpha_Intercept: ID) |    -0.72 | [-0.86, -0.46] |   100% | 1.000 | 3287.00

plot(conditional_effects(gomp.to, ndraws = 100, spaghetti = T))

pp_check(gomp.to, ndraws = 100)

re = crossing(t = seq(min(df$t), 
                       max(df$t),
                       length.out=100),
               ID = unique(df$ID)) %>% 
  add_epred_draws(gomp.to, re_formula = NULL, scale = "response", ndraws = 20)

re %>% 
  ggplot(aes(y = .epred, x = t)) +
  geom_line(aes(y = .epred, x = t, group = .draw), size = .5, alpha = .5) +
  geom_point(data = df, aes(y=l, x=t, color = ID)) +
  facet_wrap(~ID, nrow = 6, ncol = 5) + 
  scale_color_viridis() +
  ylab("Mass (mg)") + 
  xlab("Time") +
  theme_bw(12) +
  theme(legend.position = "none")

This recovers the parameter values pretty well! Note that there is some moderate correlation between \(l_0\) and \(l_{inf}\) due to sampling even if the true value is 0. The model seems to overestimate slightly this value compared to the value found in the simulated data. Let’s now plot these correlations according to the model estimates

re = gomp.to %>%
  spread_draws(# Population values
    b_l0_Intercept, b_linf_Intercept, b_alpha_Intercept, 
    # Individual offsets
    r_ID__l0[ID,Intercept], r_ID__linf[ID,Intercept], r_ID__alpha[ID,Intercept],
    # Individual variances
    sd_ID__l0_Intercept, sd_ID__linf_Intercept, sd_ID__alpha_Intercept,
    sigma) %>% 
  # Individual offsets converted onto the original length scale (in micrometers)
  mutate(L0_i = (b_l0_Intercept + r_ID__l0) * Linf_mu,
         Linf_i = (b_linf_Intercept + r_ID__linf) * Linf_mu,
         alpha_i = b_alpha_Intercept + r_ID__alpha) %>% 
  # Population averge distribution
  mutate(L0_dist = rnorm(n(), b_l0_Intercept, sd_ID__l0_Intercept)*Linf_mu,
         Linf_dist = rnorm(n(), b_linf_Intercept, sd_ID__linf_Intercept)*Linf_mu,
         alpha_dist = rnorm(n(), b_alpha_Intercept, sd_ID__alpha_Intercept))
re.mean = re %>% 
  select(.chain, .iteration, .draw, L0_i, Linf_i, alpha_i) %>% 
   mean_qi(L0_i, Linf_i, alpha_i)
  # Summarize individual values into mean, lower and upper 95 % quantiles

# Plot population average (diagonal elements)
L0_dist = re %>% 
  ggplot(aes(x = L0_dist)) +
  stat_histinterval(slab_color = "gray45", 
                    outline_bars = TRUE) +
  labs(x = "", y = expression(L[0])) +
  theme_bw(12) +
  theme(aspect.ratio=1)
Linf_dist = re %>% 
  ggplot(aes(x = Linf_dist)) +
  stat_histinterval(slab_color = "gray45", 
                    outline_bars = TRUE) +
  labs(x = "", y = "") +
  theme_bw(12) +
  theme(aspect.ratio=1)
alpha_dist = re %>% 
  ggplot(aes(x = alpha_dist)) +
  stat_histinterval(slab_color = "gray45", 
                    outline_bars = TRUE) +
  labs(x = expression(alpha), y = "") +
  theme_bw(12) +
  theme(aspect.ratio=1)

# Plot individual average with CI (lower diagonal elements)
corr1 = re.mean %>% 
  ggplot(aes(x = L0_i, y = Linf_i)) +
  geom_errorbarh(aes(xmin = L0_i.lower, xmax = L0_i.upper)) +
  geom_errorbar(aes(ymin = Linf_i.lower, ymax = Linf_i.upper)) +
  geom_point(alpha = .8, size = 3) +
  labs(x = "", y = expression(L[inf])) +
  theme_bw(12) +
  theme(aspect.ratio=1)
corr2 = re.mean %>% 
  ggplot(aes(x = L0_i, y = alpha_i)) +
  geom_errorbarh(aes(xmin = L0_i.lower, xmax = L0_i.upper)) +
  geom_errorbar(aes(ymin = alpha_i.lower, ymax = alpha_i.upper)) +
  geom_point(alpha = .8, size = 3) +
  labs(x = expression(L[0]), y = expression(alpha)) +
  theme_bw(12) +
  theme(aspect.ratio=1)
corr3 = re.mean %>% 
  ggplot(aes(x = Linf_i, y = alpha_i)) +
  geom_errorbarh(aes(xmin = Linf_i.lower, xmax = Linf_i.upper)) +
  geom_errorbar(aes(ymin = alpha_i.lower, ymax = alpha_i.upper)) +
  geom_point(alpha = .8, size = 3) +
  labs(x = expression(L[inf]), y = "") +
  theme_bw(12) +
  theme(aspect.ratio=1)

# Plot correlation estimate (upper diagonal elements)
dcorr1 = gomp.to %>% 
  spread_draws(`cor.*`, regex = TRUE) %>% 
  ggplot(aes(x = cor_ID__l0_Intercept__linf_Intercept)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linewidth = 1, color = "black", linetype = "dashed") +
  xlim(-1, 1) +
  labs(x = "", y = "") +
  theme_bw(12) +
  theme(aspect.ratio=1)
dcorr2 = gomp.to %>% 
  spread_draws(`cor.*`, regex = TRUE) %>% 
  ggplot(aes(x = cor_ID__l0_Intercept__alpha_Intercept)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linewidth = 1, color = "black", linetype = "dashed") +
  xlim(-1, 1) +
  labs(x = "", y = "") +
  theme_bw(12) +
  theme(aspect.ratio=1)
dcorr3 = gomp.to %>% 
  spread_draws(`cor.*`, regex = TRUE) %>% 
  ggplot(aes(x = cor_ID__linf_Intercept__alpha_Intercept)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linewidth = 1, color = "black", linetype = "dashed") +
  xlim(-1, 1) +
  labs(x = "", y = "") +
  theme_bw(12) +
  theme(aspect.ratio=1)

# Arrange plot into 3 x 3 grid
L0_dist + dcorr1 + dcorr2 +
  corr1 + Linf_dist + dcorr3 +
  corr2 + corr3 + alpha_dist +
  plot_layout(ncol = 3, nrow = 3, byrow = T)

Another possibility is to sample correlation estimates and plot those on a small grid with each color corresponding to the correlation strength. Solution found on Matthew Kay and Solomon Kurtz websites

levels = c("l[0]", "l[inf]", "alpha")

rho = as_draws_df(gomp.to) %>% 
  select(starts_with("cor_")) %>% 
  slice_sample(n = 60 * 60) %>% 
  bind_cols(crossing(x = 1:60, y = 1:60)) %>% 
  pivot_longer(cols = -c(x:y)) %>% 
  mutate(name = str_remove(name, "cor_ID__")) %>% 
  separate(name, into = c("col", "row"), sep = "__") %>% 
  mutate(
    col = case_when(
      col == "l0_Intercept"   ~ "l[0]",
      col == "linf_Intercept" ~ "l[inf]",
      col == "alpha_Intercept" ~ "alpha"),
    row = case_when(
      row == "l0_Intercept"  ~ "l[0]",
      row == "linf_Intercept" ~ "l[inf]",
      row == "alpha_Intercept" ~ "alpha")) %>% 
  mutate(col = factor(col, levels = levels),
         row = factor(row, levels = levels))

rho %>% 
  full_join(rename(rho, col = row, row = col),
            by = c("x", "y", "col", "row", "value")) %>%
  
  ggplot(aes(x = x, y = y, fill = value)) +
  geom_raster() +
  scale_fill_distiller(type = "div", 
                       palette = "RdBu", 
                       limits = c(-1, 1), name = expression(rho)) +
  # scale_fill_gradient2(expression(rho),
  #                      low = "#59708b", mid = "#FCF9F0", high = "#A65141", midpoint = 0,
  #                      labels = c(-1, "", 0, "", 1), limits = c(-1, 1)) +
  scale_x_continuous(NULL, breaks = NULL, expand = c(0, 0)) +
  scale_y_continuous(NULL, breaks = NULL, expand = c(0, 0)) +
  theme(strip.text = element_text(size = 12)) +
  facet_grid(row ~ col, labeller = label_parsed, switch = "y")

In conclusion

Mathematical description of the statistical models

Formally, our statistical model can be described with the following set of equations

Non-scaled model

$$ \[\begin{aligned} \large \text{Model Likelihood}\\ \text{Average growth rate}\\ L_{i,t} &\sim N(\mu_{i,t}, \text{ } \sigma)\\ \mu_{i,t} &= L_{inf_{i}}\times e^{ln \left( \frac{L_{0_{i}}}{L_{inf_{i}}} \right) \times e^{-\alpha_i t} } \\ \text{Among-individual covariance}\\ \begin{bmatrix} \mu_{L_{0_{i}}}\\ \mu_{L_{inf_{i}}}\\ \mu_{\alpha_i}\\ \end{bmatrix} &\sim MVN(\begin{bmatrix} \mu_{L_{0}}\\ \mu_{L_{inf}}\\ \mu_{\alpha}\\ \end{bmatrix}, \Sigma)\\ \Sigma &= \begin{bmatrix} \sigma_{L_{0_{i}}} & 0 & 0\\ 0 & \sigma_{L_{inf_{i}}} & 0 \\ 0 & 0 & \sigma_{\alpha_i} \end{bmatrix} R \begin{bmatrix} \sigma_{L_{0_{i}}} & 0 & 0\\ 0 & \sigma_{L_{inf_{i}}} & 0 \\ 0 & 0 & \sigma_{\alpha_i} \end{bmatrix} \\ R &= \begin{bmatrix} 1 & \rho_{1,2} & \rho_{1,3}\\ \rho_{1,2} & 1 & \rho_{2,3} \\ \rho_{1,3} & \rho_{2,3} & 1 \end{bmatrix} \\ \large \text{Priors}\\ \text{Population Intercepts}\\ \mu_{L_{0}} &\sim N(200, 40) \\ \mu_{L_{inf}} &\sim N(1500, 300) \\ \mu_{\alpha} &\sim N(0.03, 0.006) \\ \text{Among-individual variances}\\ \sigma_{L_{0_{i}}} &\sim Exp(.025) \\ \sigma_{L_{inf_{i}}} &\sim Exp(.003) \\ \sigma_{\alpha_i} &\sim Exp(170) \\ \sigma &\sim Exp(1) \\ \text{Among-individual covariance}\\ R &\sim LKJ(2) \end{aligned}\]

$$

Scaled model

Compared to the non-scaled model, the major difference here is that we put a strong prior around 1 for the \(l_{inf}\) parameter.

$$ \[\begin{aligned} \large \text{Model Likelihood}\\ \text{Average growth rate}\\ l_{i,t} &\sim N(\mu_{i,t}, \text{ } \sigma)\\ \mu_{i,t} &= l_{inf_{i}}\times e^{ln \left( \frac{l_{0_{i}}}{l_{inf_{i}}} \right) \times e^{-\alpha_i t} } \\ \text{Among-individual covariance}\\ \begin{bmatrix} \mu_{l_{0_{i}}}\\ \mu_{l_{inf_{i}}}\\ \mu_{\alpha_i}\\ \end{bmatrix} &\sim MVN(\begin{bmatrix} \mu_{l_{0}}\\ \mu_{l_{inf}}\\ \mu_{\alpha}\\ \end{bmatrix}, \Sigma)\\ \Sigma &= \begin{bmatrix} \sigma_{l_{0_{i}}} & 0 & 0\\ 0 & \sigma_{l_{inf_{i}}} & 0 \\ 0 & 0 & \sigma_{\alpha_i} \end{bmatrix} R \begin{bmatrix} \sigma_{L_{0_{i}}} & 0 & 0\\ 0 & \sigma_{L_{inf_{i}}} & 0 \\ 0 & 0 & \sigma_{\alpha_i} \end{bmatrix} \\ R &= \begin{bmatrix} 1 & \rho_{1,2} & \rho_{1,3}\\ \rho_{1,2} & 1 & \rho_{2,3} \\ \rho_{1,3} & \rho_{2,3} & 1 \end{bmatrix} \\ \large \text{Priors}\\ \text{Population Intercepts}\\ \mu_{L_{0}} &\sim N(0.15, 0.03) \\ \mu_{L_{inf}} &\sim N(1, 0.3) \\ \mu_{\alpha} &\sim N(0.03, 0.006) \\ \text{Among-individual variances}\\ \sigma_{L_{0_{i}}} &\sim Exp(34) \\ \sigma_{L_{inf_{i}}} &\sim Exp(5) \\ \sigma_{\alpha_i} &\sim Exp(170) \\ \sigma &\sim Exp(1) \\ \text{Among-individual covariance}\\ R &\sim LKJ(2) \end{aligned}\]

$$

Session info

sessionInfo()

R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=French_France.utf8  LC_CTYPE=French_France.utf8   
[3] LC_MONETARY=French_France.utf8 LC_NUMERIC=C                  
[5] LC_TIME=French_France.utf8    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] see_0.8.0              report_0.5.7           parameters_0.21.1     
 [4] performance_0.10.5     modelbased_0.8.6       insight_0.19.6        
 [7] effectsize_0.8.6       datawizard_0.8.0       correlation_0.8.4     
[10] bayestestR_0.13.1      easystats_0.6.0        viridis_0.6.4         
[13] viridisLite_0.4.2      marginaleffects_0.15.0 brms_2.20.1           
[16] Rcpp_1.0.11            tidybayes_3.0.6        patchwork_1.1.3       
[19] lubridate_1.9.2        forcats_1.0.0          stringr_1.5.0         
[22] dplyr_1.1.3            purrr_1.0.2            readr_2.1.4           
[25] tidyr_1.3.0            tibble_3.2.1           ggplot2_3.4.4         
[28] tidyverse_2.0.0       

loaded via a namespace (and not attached):
  [1] colorspace_2.1-0     ellipsis_0.3.2       estimability_1.4.1  
  [4] markdown_1.8         QuickJSR_1.0.6       base64enc_0.1-3     
  [7] rstudioapi_0.15.0    farver_2.1.1         rstan_2.26.23       
 [10] svUnit_1.0.6         DT_0.29              fansi_1.0.5         
 [13] mvtnorm_1.2-3        bridgesampling_1.1-2 codetools_0.2-18    
 [16] knitr_1.44           shinythemes_1.2.0    bayesplot_1.10.0    
 [19] jsonlite_1.8.8       ggdist_3.3.0         shiny_1.7.5         
 [22] compiler_4.2.2       emmeans_1.8.8        backports_1.4.1     
 [25] Matrix_1.6-1         fastmap_1.1.1        cli_3.6.1           
 [28] later_1.3.1          htmltools_0.5.6.1    prettyunits_1.2.0   
 [31] tools_4.2.2          igraph_1.5.1         coda_0.19-4         
 [34] gtable_0.3.4         glue_1.6.2           reshape2_1.4.4      
 [37] posterior_1.4.1      vctrs_0.6.5          nlme_3.1-160        
 [40] crosstalk_1.2.0      tensorA_0.36.2       xfun_0.40           
 [43] ps_1.7.5             timechange_0.2.0     mime_0.12           
 [46] miniUI_0.1.1.1       lifecycle_1.0.4      gtools_3.9.4        
 [49] MASS_7.3-58.1        zoo_1.8-12           scales_1.2.1        
 [52] colourpicker_1.3.0   hms_1.1.3            promises_1.2.1      
 [55] Brobdingnag_1.2-9    parallel_4.2.2       inline_0.3.19       
 [58] RColorBrewer_1.1-3   shinystan_2.6.0      yaml_2.3.7          
 [61] gridExtra_2.3        loo_2.6.0            StanHeaders_2.26.28 
 [64] reshape_0.8.9        stringi_1.7.12       dygraphs_1.1.1.6    
 [67] checkmate_2.2.0      pkgbuild_1.4.2       cmdstanr_0.5.3      
 [70] rlang_1.1.2          pkgconfig_2.0.3      matrixStats_1.0.0   
 [73] distributional_0.3.2 evaluate_0.22        lattice_0.20-45     
 [76] labeling_0.4.3       rstantools_2.3.1.1   htmlwidgets_1.6.2   
 [79] tidyselect_1.2.0     processx_3.8.2       GGally_2.1.2        
 [82] plyr_1.8.8           magrittr_2.0.3       R6_2.5.1            
 [85] generics_0.1.3       pillar_1.9.0         withr_2.5.2         
 [88] xts_0.13.1           abind_1.4-5          crayon_1.5.2        
 [91] arrayhelpers_1.1-0   utf8_1.2.4           tzdb_0.4.0          
 [94] rmarkdown_2.25       grid_4.2.2           data.table_1.14.8   
 [97] callr_3.7.3          threejs_0.3.3        digest_0.6.33       
[100] xtable_1.8-4         httpuv_1.6.11        RcppParallel_5.1.7  
[103] stats4_4.2.2         munsell_0.5.0        shinyjs_2.1.0